#=====================================================================
# 2020/07/27
# Geothermal prediction dataviz
# Fumiya Uchikoshi, uchikoshi@princeton.edu
#=====================================================================
rm(list=ls())
######################################################################
# Loading packages
######################################################################
library(gdata)
library(tidyverse) 
library(dplyr)
library(ggthemes)
library(knitr)
library(readxl)

######################################################################
# Assign colors
######################################################################
cbp2 <- c("#999999", "#E69F00", "#56B4E9", "#009E73",
          "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

######################################################################
# Onsen between
######################################################################
onsen <- sort(c(rep(c(0,15,30,45),2))) 
onsen <-  as.data.frame(onsen)
df <- read_csv("../Results/Prediction_OnsenBtw.csv",skip=1) %>% 
   filter(X1 == "b" | X1 == "se") %>% 
   bind_cols(onsen) %>% 
   pivot_longer(ends_with("yearx"), names_to = "year", values_to = "value") %>% 
   pivot_wider(names_from = c("X1"), 
               values_from = value) %>% 
   mutate(year = as.numeric(substr(year,1,1)),
          year = year*5+1975) %>% 
   mutate(onsen = paste0(onsen," ryokan")) %>% 
   mutate(yearx = year + 4,
          year = paste0(year,"-",yearx)) 

df %>% dplyr::select(1:3) %>% 
   write.csv("../Results/Figure2.csv")

df %>%
   mutate(year=factor(year,levels=unique(df$year)),
              onsen=factor(onsen,levels=unique(df$onsen))) %>% 
   ggplot(aes(x = year, y = b, color = onsen,group=onsen)) + 
   geom_point(aes(shape=onsen))+geom_line(aes(linetype=onsen)) + 
   ylab("Probability of geothermal drilling exploration (0-1 scale)") + xlab("Year") +
   theme_few()+scale_color_manual(values=cbp2)+theme(legend.title=element_blank(), legend.position = "bottom",
                                                     axis.text.x = element_text(angle = 30, hjust = 1))
ggsave(height=6,width=9,dpi=200, filename="../Results/Figure2.pdf",  family = "Helvetica")

######################################################################
# Farming within
######################################################################
year <- sort(c(rep(c(0:7),2))) 
year <-  as.data.frame(year)
df <- read_csv("../Results/Prediction_FarmingWth.csv",skip=1) %>% 
   filter(X1 == "b" | X1 == "se") %>% 
   bind_cols(year) %>% 
   pivot_longer(ends_with("_at"), names_to = "farming", values_to = "value") %>% 
   pivot_wider(names_from = c("X1"), 
               values_from = value) %>% 
   mutate(year = year*5+1975,
          farming = as.numeric(substr(farming,1,1)),
          farming = case_when(
             farming==1 ~ "Town C: Start at 28% above its average farming employment over time, decline 8% per period",
             farming==2 ~ "Town B: Start at 14% above its average farming employment over time, decline 4% per period",
             farming==3 ~ "Town A: Start at 0% above its average farming employment over time, decline 0% per period")
          ) %>% 
   arrange(farming)%>% 
   mutate(yearx = year + 4,
          year = paste0(year,"-",yearx)) 

df %>% dplyr::select(1:3) %>% 
   write.csv("../Results/Figure3.csv")

df %>% mutate(year=factor(year,levels=unique(df$year)),
              farming=factor(farming,levels=unique(df$farming))) %>% 
   ggplot(aes(x = year, y = b, color = farming,group=farming)) + 
   geom_point(aes(shape=farming))+geom_line(aes(linetype=farming)) + 
   ylab("Probability of geothermal drilling exploration (0-1 scale)") + xlab("Year") +
   theme_few()+scale_color_manual(values=cbp2)+theme(legend.title=element_blank(), legend.position = "bottom",
                                                     axis.text.x = element_text(angle = 30, hjust = 1))+
   ylim(0,.04)+
   guides(col = guide_legend(ncol = 1))
ggsave(height=6,width=9,dpi=200, filename="../Results/Figure3.pdf",  family = "Helvetica")
